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I. INTRODUCTION 


The basic factors controlling the strength of materials are found 
at the atomic scale of material structures. The availability of this type 
of structural detail is essential to the interpretation of available 
experimental data. It lies at the heart of our ability to predict and 
modify the strength of materials to suit our technological needs. 
Unfortunately, the direct experimental observation of the atomic scale 
phenomena determining the strength of metals has proven impossible; other 
approaches are needed in order to obtain the atomic scale picture. We 
at Battelle have developed a comprehensive mathematical model by which the 
collective behavior of a very large number of atoms within a metal or alloy 
can be accurately simulated^^\ In particular, the manner in which the 
atomic interactions relate to dislocation motion and crack extension to 
determine the strength of a metal can now be determined and dissected by 
this method. Our current effort is in understanding the factors determin- 
ing resistance to crack extension in iron metal subject to a stress, and 
how the presence of hydrogen causes crack extension to proceed at 
abnormally low stresses- 
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II. TECHNICAL PROGRESS 


A<p Overview of the Areas of Current Progress 

■When this project started its main emphasis was on th*’. calcula- 
tion of the interatomic forces. This is one of the basic ingredients in 
the simulation of stress crack extension of metals in the presence of 
hydrogen. However, in the process of exercising the crack simulation 
model, even with empirical interatomic forces, it became clear that some- 
thing was seriously wrong. The effort in the calculation of the inter- 
atomic forces was consequently reduced, and attention was focused on 
Identifying the source of the problem in the crack simulation model. The 
material in this report thus reflects the spread of the effort to both 
complementary aspects of the project. The material is presented in the 
following order. 

e Brief sketch of the problems encountered with GehXen’s 
crack simulation model 

« Identification of the source of error in the crack 
simulation model 

o Hew improved formulation of the crack simulation 
model 

o Electronic structure calculation of small iron atom 

clusters, their stability as a function of configuration, 
and the effect of adding a hydrogen atom. 

B . Sketch of the Problems 


The basic crack simulation model with which we have been working 
is that pioneered by P, C, Gehlen at the Battelle Columbus Laboratories. 

In order to place Che magnitude of the difficulties with the model in 
perspective it is necessary to present first some background information. 
The physical problem being studied is the resistance to crack extension 
which arises from the lattice structure in metals. This resistance to 
crack extension is shown in Figure 1 in the form of the energy barrier 
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The and AEj curves cross at what is called the Griffith stress intensity'' * . 
At this point the rate of continuum strain energy release equals Che rate 
of surface energy j , creation 

for 

where E is Young’s -mcxiulus. Although from an elasticity-theory point of 
view cracks become unstable for K^Kg, these do not propagate readily 
because of the ’’lattice trapping” illustrated in Figure 1 as first pointed 
out by Hsieh and Thomson^ - Only for K “ K_ does all resistance to crack 
extension dissappear as the barrier goes to zero. 

Contrary Co what one expected, the application of the crack 
simulation model to this problem predicted that the barrier to crack 
extension would not disappear. It led to results of the type shown in 
Figure 3 . 





Fi&u^e 


This unphvsical result was attributed to the use of a linear elastic field 
in the description of the continuum region. The latter is the region within 
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which the discrete crack core region is embedded. Subsequentl;^ , a consider- 
able amount of effort was expended in replacing the linear elastic field 
by a general non-linear field which was based on the use of Green's functions 
to fold residual forces into modified equilibrium displacement fields^^^. 
Unfortunately, when the corresponding energy barriers as a function of 
stress intensity were calculated by this considerably more sophisticated 
method, they made even less physical sense than those obtained previously. 

Furthermore, disturbing convergence to different relative 
mluima were found depending on the procedure used to locate these 
minima and convergence to the saddle points was plagued with difficulties. 

The magnitude, persistence, and uncertain origin of these problems, in 
spite of the constant and well reasoned efforts to overcome them, persuaded 
P. C. Gehlen, the originator of the simulation model, that the approach 
lacked validity. He was discouraged from further pursuing the source of 
the difficulties. 

C. Identification and Analysis of the Sources of Error 


We have given Gehlen* s procedures a detailed examination. This 
effort has brought to light some fundamental assumptions in his model which 
are found to be incorrect. The deficiencies in Gehlen 's model, in 
essence, fall into three categories although they are interrelated. 

o Dis registry between the continuum and discrete region 
crack fields 

o Neglect of the contribution of the continuum strain 
energy release 

« Inconsistent definition of potential energy and 
gradients thereof 

We shall discuss each of these points next. 

1. Disregistry Between the Continuum 
and Discrete Region Crack Fields 

We found that in the simple rigid boundary model, apart from 
restricting the elastic crack field to be linear, the origin of the field 
was always held fixed at an arbitrary position. Thus, even when the 
discrete region would describe an advanced crack, e.g. saddle point or an 
adjacent stable crack, tne field describing the continuum would still be 
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fixed in its arbitrary original position. We refer to this as the 
disregistr y effect arising from the incorrect relative alignment of the 
fields in both discrete and continuum regions. Ihis is depicted 
schematically in Figure 4. 



F\EI,D3 m 
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It was found that when the disregistry was minimized by advancing the 
elastic fields together with the discrete region, the relative stability 
of two adjacent equilibrium cracks changed by amounts of the same 
magnicude as those of the energy barriers we were looking for. Thus., 
the disregistry effect was found to be an important source of the errors 
in the results shown in Figure 3. The introduction of the flexible 
boundary model in principle solves this problem. It is designed to fully 
relax the continuum region, iiO- to both introduce non-linear effects and 
to minimize the disregistry effect. In practice the relaxation of this 
elastic field w .is not direct , but entered only as a response to a separate 
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full relaxation of the atoms in the discrete region. The latter region 
was relaxed hy the “ kinetic energy quenching” method^^^ . There remains 
In practice the possibility that, as the kinetic energy is quenched in 
the discrete region, the rate at which information about disequilibrium 
la transferred to the continuum region becomes exceedingly small. This 
would result in only a partly minimized disregistry effect. This 
however, was probably only a minor error in the flexible boundary model. 

The major reason why this improved model failed was an imcomplete account- 
ing of the terms in the total energy of the system. This is discussed 
in the next suhsection. 

Before going on it is worth noting aspects of the saddle point 
calculations which were affected by the disregistry. In the rigid boundary 
model, the same type of disregistry error was introduced as in searching 
for equilibrium minimum positions. In the flexible boundary case, the 
continuum atoms were left positioned at the maximum of the energy along a 
straight line interpolating from their position in the t\JO equilibrium 
cracks on either side of the saddle point. Although the continuum was 
not further relaxed in searching for the saddle point, this is probably 
by itself only a minor effect. The major error in the saddle point 
calculation arises from the calculation of forces on the atoms by a 
formula which is inconsistent v.ith the definition of the potential energy 
of the model. This will be discussed in the third subsection. 

2, Neglect of the Energy Exchange 
With the Continuum 

When the elastic field of the continuum region advances in 
order to minimize the disregistry effect with the discrete region, then, 
from elasticity theory, one can calculate the strain energy released from 
the continuum which flows into the crack tip region as depected 


schematically in Figure 5. 



Figure 5. 


This strain energy release 


( 2 ) 


is given by 


E 


where K is the stress intensity, E is Young’s modulus, Ab is the advance of 
the crack elastic field, and "a" stands for the thickness of the slices 
in the model. This energy change of the continuum has to be combined with 
the energy change of the discrete region. These two energy changes are 
shown in Figure 6. 
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The energy change of the discrete region is typical of calculations where 
the disregistry effect has been minimized by advancing the continuum together 
with the discrete region crack. At a repeat distance Ab one finds that 
this energy has increased (within the accuracy of the model) by the surface 
energy created, i.e. 2'yaAb , where "a" stands for the width of the main slice 
in the model- The combined energy change is shown in Figure 7 for the 
case when crack growth is favored. 
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In Gehlen's rigid boundary model the disregistry effect originating from 
not allowing the continuum field to advance resulted in an artifactually 
zero continuum strain energy release. The combined consequence of these 
restrictions is the underestimation of the basic driving force behind the 
disappearance of the effective energy barrier to crack extension, as is 
seen from Figure 7. This explains why the results of Gehlen’s rigid 
boundary calculations shown in Figure 3 predicted that the barrier to crack 
extension would not disappear. In Gehlen’s flexible boundary model the 
disregistry effect was removed by allowing the field to readjust. The 
calculation of the strain energy changes associated with these continuum 
readjustments should provide the required driving force lacking previously. 
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We £ound, however, that these continuum energy changes were unfortunately- 
not Included in Gehlen's calculations. The omission of this major energy 
contribution in both the crack calculations of m/.nima and caddie points 
is probably the principal reason for the unphysical nature of the results 
obtained from the flexible boundary model. 


3. Inconsistent Definition of Energy and 
Gradients of the Energy 

The forces on each of the 'atoms are the central ingredient in 
the search and location of the equilibrium crack minima and saddle points. 
The definition of the force on the k-th atom used in Gehlen's model is* 



— V u 



This is also the definition which makes most physical sense. The potential 
energy in Gehlen's model is defined in turn as 





tc 


; 0 . 
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containing half of the interaction energy with the continuum. The objective 
definition of the force on an atom is the negative gradient of the 
potential energy. Thus, the definition of the model potential energy fixes 
the definition of the force. In this case we find 



It is clear that this is not the definition of the force used by Gehlen. 
It is not the natural definition since it weights the continuum atoms' 
contributions differently than the discrete region atoms, yet it is the 
definition which is consistent with the properties of the model as 
contained in the potential energy. Gehlen's definition of force would 
be consistent with a potential energy of the form 


= f L S u(IR,-Rtl) ^ S E, -iilVUi) 


keO 




See Section F, page 22, for definitions of the terms used here, 
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The use of Gehlen's definition of force, Ik , in the relaxation techniques 

( 2 ) 

leads to the minima and saddle points of E .If this were also the 
potential energy used to calculate the associated energy differences, then 
there would be no inconsistency error. In Gehlen’s model, however, 
was used in calculating the energy difference. It was found that the 
equilibrium configurations of and E^^^ could differ significantly. 

This is thus another major source of error in the simulation model. It 
permeates all the applications and versions of the simulation model. In 
a manner, the inconsistency problem’ arises because not all necessary 
energy contributions of the model were built into the potential energy. 
This is just another aspect of the neglect of the energy terms from which 
we obtain the continuum strain energy release. 

4, Acknowledgements 

The analysis and identification of the errors discussed in this 
section are the result of extended discussions by the author with 
Drs. P. C. Gehlen, G. T. Hahn, R. G. Hoagland, and A. J. Markworth of 
the Battelle Columbus Laboratories Metal Science Section. 


D. Improved Formulation of the Crack Simulation Model 


Our analysis of Gehlen’s simulation model shows that, apart from 
the three previously described conceptual errors, the remaining aspects of 
his model retain their validity. We have built upon Gehlen's approach in 
formulating an improved simulation model which incorporates the solution 
to the conceptual errors in the previous model. The details of the new 
approach are described in the Appendix to this report. The main points 
are 

« The total energy of the system has been formulated precisely 
with emphasis on including the contribution from atoms in all 
regions o<^ the model. The result is 


E - E, + S ! + + ^2 




The term E^ is the energy of the discrete region atoms 
defined as ( T is a translation vector; see Appendix) 

r*- U 
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e The total energy of the system has been formulated precisely 
with emphasis on including the contribution from atoms in all 
regions of the model. The result is 


E = E, + 


il£ p 0 


The term is the energy of the discrete region atoms 
defined as (T. is a translation vector; see Appendix)- 

fV-' iti 



'h ^0 

f, [udar^tH + Z 

ktb '^=<4 

The second term is a work energy term giving the energy change 
o£ continuum region atoms relative to a reference strained 

O® r- 

State ^ with energy . This energy is the source of the 

continuum strain energy release, , discussed in 

subsection C.l . The gradient of in the work energy is 
defined as 



e The formulation of the complete energy of the system allows 
us to unify the treatment of the discrete and continuum 
regions in the search for equilibrium. The search for 
equilibrium becomes a problem of finding the stationary 
points of this total energy \-7ith respect to all the degrees 
of freedom in the model. The problem of minimizing the 
disregistry effect become.^ simply part of the process of 
seeking for the lowest total potential eniirgy, 

o The application of the "con'jugate gradient" technique of 

( 5 ) 

Fletcher and Reeves is discussed as a natural approach in 
light of the function minimization mold into which we have 
now cast the simulation model. 

o The previously used relaxation technique of "kinetic energy 
quenching" is extended to the direct relaxation of the 
parameters of the elastic continuum field. 

e The forces on the degrees of freedom of the model are 
found directly from the gradients of the total potential 
energy definition. 

It is worth noting that in the essential aspects this improved 

( 6 ) 

formulation brings Gehlen’s basic approach closer to that of Sinclair' 

One remaining difference lies in the now vivid perspective we have on the 
singular role the present modifications of Gehlen’s model have in making 
this a valid approacht 
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£• Electronic StructuTe Calculations of Clusters of Iron 
Atppus anii The Effect of Inserting Hydrogen 

Our goal is to obtain detailed information about the Interaction 
of a hydrogen atom with an environment of iron atoms. This is the basic 
Ingredient to the crack simulation of hydrogeu'-enhanced stress cracking of 
iron. Our approach is to carry out a series of cluster calculations with a 
hydrogen atom placed variously within iron atom clusters. An example 
configuration is that of hydrogen at an octahedral site in bcc Iron as 
shown in Figure 8. 



These calculations will give direct information about the electron density 
redistribution in the presence of hydrogen. Moreover, from the calculated 
elements of the electronic energy surface, we expect to construct a rapidly 
convergent many— body decomposition of the form 
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which contains the important geometry dependence of the total energy 
surface. This is the direct ingredient to the crack simulation model. 

The electronic structure calculations are being carried out in 
collaboration with Dr. Richard Jaffe of the NASA Ames Research Center. 

They en^loy a novel technique for folding the effects of core electrons 
Into the form of effective pctentlals which makes the cluster calculations 
tractable by the standard methods. The latter techniques are Implemented 
using the very efficient GAUSSIAN 70 code. Finally, the magnitude of 
these calculations requires the special con^utatlonal facilities available 
at NASA Ames. 

We started by first checking the use of the UHF Code in dealing 
with the high spin cases posed by iron atom. The ground state, can be 
represented in terms of real orbitals in a single determinant wavefunction 





1 * 

r^r 

L 







which contains four unpaired orbitals. The energy for this wavefunction, 
and various excited states, was studied as a function of the contraction 
in the primitive gaussian basis consisting of a 3s, 3p, and 5d functions. 
The results are shown in Figure 9. These results show that there is little 
loss in excitation energies in contracting the (335) set to [312], This 
alone seems, however, to be a stricter criterion of the quality of a 
contraction than is necessary in calculating the ground state atom-atom 
interaction. We find, as shown in Figure 10, that the Fe-Fe potential 
energy curve is essentially Identical for a [311] set as for a [312] set. 
Even the calculations in a minimal [llll.set yield a reasonable Fe-Fe 
potential curve in this ease. 
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v+ 

The Icwest elecCronlc state was found to be a 2^p state. In 
this state two Iron (L electrons remained singlet paired, and the remaining 
ones all align in the same direction. This high spin electron configura- 
tion Is depicted schematically as 


Hiw SPIM ; Fe ['♦iM 3d''llH'l1l)] + Fe[di[lF)3i‘’(lH'm)] 

s V ^ 

t t 

The doubly occupied dt orbital was found to be a d^ orbital where that axis 
^Incldes with the intemuclear axis. On the basis of electron repulsion 
arguments one would have predicted this orbital to be d^. The low spin 
state was compared against the high spin state at R « 5.404 bohr. The 
electron configuration is depicted schematically as 


bw •. 


^ ^ 


+ Fe[4i(n)3d’(rv FF h)] 



The energy of this state was found to be 0.096 Hartree/atom more repulsive 


than the high spin state at the same geometry. 

Similar calculations were made for Fe^ where all Fe atoms are in 
a plane, the [311] Gontraction was used since it had been found satisfactory 
in the Fe 2 calculations. The results are shown in Figure 11. In contrast 
to the Fe 2 case, we find that in Fe^ the low spin state rather than the 
high spin state is the energetically lowest one. The intermediate spin 
state lies energetlGally between these two limits. It is Interesting to 
note that as in the Fb 2 case, the Fc^ energy as a function of a symmetric 
stretch is purely repulsive. The introduction of correlation effects will 
probably at least yield a Van der Waals minimum. The major effect in 
decreasing the repulsive energy will probably arise just frcjm changing 
the geometry from planar to tetrahedral since the number of bonds between 
iron atoms then Increases from 4 to 6. 


The most interesting results were obtained upon introducing a 
hydrogen atom into the four iron atom cluster. The results are shown in 
Figure 12. The presence of the hydrogen atom stabilizes the iron cluster 
despite .’ven the repulsive iron-iron interactions. It is surprising that 
such a large effect arises from just one additional atom. Moreover, the 
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•ffecc of the hydrogen seems to be a short range effect since the cluster, 
energy does not change until the iron atoms are about 5.6 bohr from the 
hydrogen atom. We also note that the high spin iron cluster becomes the 
lowest energy state when Interacting with a hydrogen atom with spin 

aligned in the same direction. The opposite combination of spin alignments 
produces little change over the Fe^ high spin Cluster results. Insofar as 
the electronic structure which xjnderlies these results is concerned, we 
find that there is a net charge transfer onto the hydrogen atom from the 
Iron atoms. This is shown in Figure 13 by indicating the net nximber of 
valence electrons on each atom. 
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re — 



Figure 


We see that the hydrogen atom has a 0.56 increase in electron charge. The 
process of charge transfer is actually found to be a combination of effects. 
The (j(. spin orbitals show the hydrogen atom donating 0^27 electron charges 
into iron 4p orbitals, while the spin orbitals show the hydrogen atom 
accepting 0.83 electron changes from the iron 4s orbitals. Thus, charge 
transfer is an important electronic factor in the bonding of hydrogen 
to the iron atom cluster. Preliminary results indicate that, in contrast, 
a helium atom produces almost no energetic change in the iron atom cluster. 
This Is quite reasonable given the stability of the helium atom with respect 
to donating or accepting extra electron charges. 
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Flnallyt we obtain also the energy surfacf^ ior the motion of 
hydrogen atom about a fixed configuration of four Iron atoms. The 
energy contours are shown in Figure 14, The most stable position is at 
center of the cluster* The energy increases most rapidly in the 
direction of the iron atom. It Is interesting to compare this against 
the results predicted by use of the pair potential superposition* We 
show two different empirical Fe-H pair potentials InvFigure 15 together 
i^th the Fc—Fe interaction for reference. In Figures 16 and 17 we show 
the corresponding total energy contours resulting by use of these pair 
potentials* The energy surface In Figure 16 most closely resembles the 
present results of Figure 14. However, the predicted binding energy Is 
three times too large. The energy surface in Figure 17 predicts that 
hydrogen will not be bound even at the center of the cluster* It is 
least in accord with the present results. 

F* Definition of Terms in Formulae of Subsection C.3. 
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set of all atoms in the discrete region 

set of all atoms in the continuum region 

position vector of the k-th atom 

gradient operator in the *•*.. coordinates 

interaction potential between any pair of atoms, and 
dpendent on just the distance between them^ " . 
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III. DIRECTIONS FOR FUTURE WORK 

Dislocation motion is among the most Important accepted mechanisms 
for the plastic deformation of metals in response to an applied stress. 
Gehlen found that the simulation model predicted that an edge or screw dis- 
location would not move even at abnormally high stresses . Our analysis of 
Gehlen 's simulation model finds errors which are few but which fundamentally 
af fect all results - derived from its use* The improved formulation removes 
these errors from the model. It is thus now possible to reexamine these 
basic questions about dislocation motion with the improved model. With 
little additional effort, it is also possible to examine the manner in 
which an impurity atom, such as hydrogen, ping a dislocation. Beyond 
these questions, there are basic phenomena related to dislocation inter- 
actions, such as dislocation pile-ups and annihilations, which can 
readily be treated with this model. 

The improved simulation model allows us also to resume on a firm 
basis the study of the threshold for crack propagation in iron and of the 
effect of hydrogen thereon. Although the errors in Gehlen' s model had 
been identified and verified in essence, the correction factors have 
thus far been introduced only in an a posterigri manner. The successful 
solution of the dislocation and crack motion problems, will enable us to 
consider for the first time the crack propagation problem in the presence 
of a dislocation. The latter will introduce in our simulation the 
competing mechanism of crack energy dissipation by plastic deformation 
via dislocation motion. 

Thu electronic structure calculations should explore the effect of 

additional shells of iron atoms, possibly with use of effective potentials 

to fold away the effect of all ci shell electrons bn outer shell iron atoms. 

In general, we need to develop and test rigorous systematic procedures 

for folding the external cluster effects, including approximations 

thereto , into the equations determining the electronic structure inside 

the clusters. We also should explore the introduction of correlation 

effects beyond the UHF approximation. One promising approach is the many- 

( 8 ) 

body perturbation method as used by Bartlett^ . A separate but urgent 
problem is to develop functionals for two and three body potentials which 
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lead to reliable local maay~body decompositions of the total energy 
■urfece. In this connection, it may be worthwhile to explore directly 
■eny-body decompositions of the electronic energy gradient as an alternative 
approach. 
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A. Definitions of renlons of the model 

We start by defining the different regions In the model. For this 
purpose we refer to Figure 1 first. 



The Immediate vicinity of the defect, e.g., the crack tip. is surrounded by 
the cylindrical region shown above, and which is such that all points within 
this region have coordinates that satisfy the relationships 

« i. X*” + i ij' 


This region is to be called the discrete region, U. All atoms assigned to 
this region are to have no restriction on their optimal position except for 
those dictated by the interatomic forces. 






1hl« region Is to be called the continuum region C. All atoms assigned to 
fchta region have positions dictated by elasticity theory. 

Finally, there are the two regions adjacent to C and 0. All points 
in these regions satisfy the condition 


li 1 \ 


All atoms assigned to these regions have positions fixed by translational 
synnetry* The translation vector Is written as 


0, 


-n - 0, tl , ± 2 ^ 1 , 


The atom positions are obtained by 




k 7 >1 t 


= /I + I 4 1 -4- 


The n “ 0 case gives the atoms in the main slice . The n ft 0 give each 
the atoms in the adjacent slices . This is schematically illustrated in 
Figure 3. I'V 
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B* The configurational energy 

For .the preseat discussion we will assume that the configurational 
energy Is given by the pair-wise Interactions of all the atoms, l.e., 

E - t II it(l£rSel) 

< I 

Furthermore, the configurational energy we shall be discussing is that of atoms 
In one slice, e*g. , the main slice (MS), rather than the energy of all the 
atoms. In order for the configurational energy to be the same for each slice 
ve shall include in the energy per slice half of the interaction energy with 
the adjacent slices (AS ) . 
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Upon subdividing the atoms in each slice into atoms in the D and C regions 
we obtain further 


E ' 
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This energy expression contains terms where the relevant part is embedded 
in a divergent series. These terms require separate analysis to follow 

later. 

The explicitly convergent terms correspond to the energy of the 
discrete region and its interaction with the continuum. One can readily 
show that 


-o'j 


i l l r, U(|R,-Rj-Tj| 
Uii 11=^ ’’ 




tI I I 

•v.^f 


u 



ttslng the ^relation 


T 

-w-n. 



Uelng this result one may write the explicitly convergent terms as 


+ 



■ tv A. < - • J J 


the type of interaction appearing in the first summation Is represented 
schematically In Figure 4. 



FIGURE 4. 

Similarly, the second summation contains interactions shown schematically 
in Figure 5. 



FIGURE 5. 
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Note that the finite range of the pair potential allows one to 
truncate the summation over adjacent slices to Just the first pair of 
adjacent slices. The summation over continuum atoms can also, for the 
same reason, be truncated to the finite number of continuum atoms which 
are within the range of the pair potentials of the discrete region atoms. 

He shall refer to this subset of the C region atoms, as the C* region atoms. 

The divergent terms in the energy correspond to the energy of the 
continuum region. These terms are 






The type of interaction appearing in these sums is represented schematically 
In Figure 6. 



FI6UBE 6. 


This energy is Infinite in the stress-crack, problem since it corresponds to 
the total strain energy of the body, and this energy diverges as the size of 
the body increases indefinitely. However, the change in the strain energy as 
the elastic crack field is modified is finite. Therefore, we consider the 
energy of the continuum configuration, ^ i’ , relative to the ejiergy of a 
reference strained continuum configuration, r/ t t- , which is Itself also 
infinite, but relative to which the energy change is finite 





f 4fl i- 

X I' 


The change in E„ is a path integral of the directional derivative along a 
path Joining the two configurations r. and > . This Integral gives the 
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vork energy involved In moving the system from configuration ^ to con- 
figuration % . 

This path integral is independent of the choice of path. We 
choose a conveniently simple path, namely, a straight line joining the two 



The locus of all points along this path is given by the vector 


1 == li* * 

Along this path we have 


» /V j 


i X i 
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Using these results in the path integral we obtain 
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nius, we have now resolved the change In the strain energy of the continuum 
Into the contribution from each of the atoms In the continuum region. The 
■omnatlon over the continuum region atoms truncates to just those continuum 
atoms which are within the range of the pair potentials of the discrete 
region atoms; l.e., the atoms in the C* regions. The contribution of the 
continuum atoms not In C* vanishes because for these atoms 


- V 

and 

V = 0 


from the equilibrium condition for all points in the continuum displaced 
according to elasticity theory. 

In summary, the configurational energy per slice Is now specified 
precisely as 



'T'f'i 


We shall not at present assume the equilibrium relation • 

Approximate formulae for the work energy term may be derived by 
use of the trapezoidal rule and Simpson's rule for the calculation of Integrals. 
Define 


f 



% 
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Then applying che trapezoidal rule one finds 

I --t[F^(R) + F,(n] 


and applying Simpson's rule one finds 


I 

I ViEjUliX ■ ■■f F|.( + i J') + ] 


C« Validity of the work energy formula 

In truncating the work energy in E, 


T. (£rL;).i \E,'.:oA - E 


b 0 

we have assumed that 




for all atoms In C which are not in G*. This is justified to the extent 
that the atom positions are given by elasticity theory, and that elasticity 
theory is valid in this region. We shall ass\jme the latter since this can 
always be made the case by enlarging the size of region D. Insofar as the 
atom positions are concerned, it is basic to the model that the continuum 

e* 'I 

configuration, r 1: , Is given by elasticity theory. Moreover, the reference 
continuum configuration, , can always be chosen to satisfy elasticity 

theory. Thus, the Integrand 


vanishes for A.= 0 and X = 1 . It Is not apparent that it also vanishes 
for the values of \ occurring in between, 0 \ a < 1 » We shall show 
that In the case that the crack tip moves in a straight line joining the 

and I limits, then, within the approximations of linear elasticity, 


the Intermediate configurations « 0 ^ X O t fact equilibrium con- 

fi 4 |urations. 

We start by specifying the end points as equilibrium configurations 
by writing 

'll' * t ni(r(f-t’] 

“V nJ f\, \ I 

^ere are the perfect crystal positions of the atoms, £ is the posi- 

tion of the crack tip, and stands for the displacement field given by 
linear elasticity. All intermediate equilibrium configurations are obtained 
from 


i - C + Kiajf-'-t' 


note that 


w( t ] = U/( r"- 




where 


lit - t-t 


/w r\f 


Similarly, we also have 
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where 


At' - t - 1’ 


Therefore we can write both 


= R.’ + lA vur( 'il'’'- 1’ ,At' 4 


and 


« = r' + k vu;( f - 1' ). At’ 4 „ 


Next, we choose the crack tip to arrive along a line joining the endpoint 
crack-tip positions 


t = t' + \ ( t' - t' '■ 

Kj rV \ ^ 


fsj IV 


0 Ta < ^ 


Using this iii the equations for f , and expanding 7 tC, t,' ) about 

'R^tM one finds 
1 rj 'V I 


ft - R“ 4 XK VUJlf-t'l.lt'-t': 4 ... 

r\f <\y i ' 


and 


ft - a' + (x-on7uif-t’).it'-t”) + ... 


Neglecting all terms which are second and higher order changes In , as 
is consistent with linear elasticity, and subtracting both equations, after 
first dividing each by \ and respectively, we obtain 
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or 


ft = ft" X ( ft' - ft’ 1 

^ ^ I A/ J 


This shows thatt within the assumptions of linear elasticity, the configurations 
generated by passing a straight line between two equilibrium configurations 
are themselves equilibrium conflgruations associated with a crack tip which 
is also moving on a straight line. 

D. A new approach to the search for the equiljbriiiin configuration 

One of the important lessons from our past work is that the continuum 
elastic field, in which the discrete region is embedded, cannot be left posi- 
tioned arbitrarily. The effect of variously positioning the elastic field, 
followed in each case by a full relaxation of the discrete region, leads to 
energy changes of the very same magnitude as the energy differences we are 
interested in. It is clear that the position of the cgatinuum field has to 
he carefully adjusted in direct concert with the relaxations occurring in the 
discrete region. The aim is to minimize the disreglstry between the fields 
in the two regions. 

The configurational energy (potential energy) of the model, whicdi 
was discussed in section B, depends directly on the position of the con- 
tinuum field through the positions of the atoms in this region. It is there- 
fore possible to attain minimal disregis try simply as part of the process of 
seeking for the lowest total potential energ y in all variables, i.e., includ:Jng 
the position of the continuum field with the other degrees of freedom relative 
to which the system is being relaxed. This approach affords a unified treatment 
of the approach to equilibrium of the discrete and continuum regions. This 
unity of approach extends to the introduction of non-linear effects in the 
elastic field. In this more general ease, the parameters adjusting the char- 
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ACter and admixture of non-linear elastic field terns are simply to be added 

(In the list of variables to be relaxed) to the position of the overall field 

and the positions of the atoms in the discrete regiot: . 

In the pioneering ulmulatlan work, of F. C. Gehlen' , the discrete 

and continuum regions were cast in a mold which naturally pointed to separate 

and distinct procedures for relaxing each region* Thus, the discrete region 

( 2 ) 

was relaxed using the technique of * *klnetlc energy quenching" by Larsen' , 
while the continuum was relaxed (In the flexible boundary version) by 
generating, via Green’s functions, the displacement fields which drove 
all the current non-zero forces to zero. These two processes were used 
sequentially until the system was found In equilibrium. 

(3) 

In the present approach, as also in the simulation work of Sinclair , 

all distinct properties of each region are built in at the outset in the con- 
figurational energy functional. Thereafter, the search for equilibrium be- 
comes a problem of minimization of this functional with little remaining 
distinction between the discrete and continuum regions. One approach which 
seems practical for a function minimization involving as many independent 

variables as found in the crack simulation (300 - 1000) is the method of 

(4) 

” conjugate gradients' * of Fletcher and Reeves . This method requires in 
essence only the calculation of the gradients of the potential energy in 
order to lead to a minimum. It involves a sequence of search direction 
vectors (constructed from the gradient informatioa) and one-dimensional 
function minimizations in these directions. Define the gradient vector as 

& - VE 

where E is the configurational energy functional defined in section B, and 
the gradient vector has a component for each linearly independent variable 
in E. The Fletcher and Reeves algorithm is 
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This function talnimlzation method also has the advantage of requiring but 

relatively minor modification to serve in finding saddle points. It then 

C5) 

becomes the method of Sinclair and Fletcher . A closely related function 
minimization algorithm is that of Fletcher and Powell^^\ Although this 
latter algorithm may be generally more stable and rapidly convergent there 
is a price. It requires at each step the additional construction and mani- 
pulation of a full (symmetric positive definite) matrix which at convergence 
becomes the Inverse of the Hessian matrix. Moreover, one requires 
additional storage spaces, where M is the number of independent variables. 
This limits the application of the Fletcher and Powell method to cases for 
which . 
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A method which has been successfully used in the past to find 
•qullibrlum configurations is the "kinetic energy quenching " method of 
Leraen^ . This method is based on the integration of the classical 
equations of motion as a convenient device to enable the atoms to probe 
the configurational energy surface. This Is coupled with periodic re- 
movals of the accumulated kinetic energy until the system Is at a potential 
energy minimum. While the classlcal-equations-of-motions method is 
straight forward when applied to the motion of the atoms in the discrete 
xeglon» such is not the case when, as in the present approach, we also 
wish to relax by the same procedure the parameters of the continuum elastic 
field. We show next that Indeed it is possible to extend this method to 
deal directly with the elastic field. 

The kinetic energy of the system is the kinetic energy of the 
atoms in the discrete and continuum regions respectively 


K = 1: 


. 'V r 
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We constrain all points of the continuum region to be displaced according 
to an elastic field specified by a set of m parameters, /; . . i 

1* €t • I 





We define the mass tensor 








The kinetic energy may now be written as 
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vh«re the kinetic energy contribution from the continuum atoms is now 
expressed In terms of the mass tensor and the time rate of change of the 
parameters of the elastic field. Next we define the momentum associated 
with each of the degrees of freedom. 



«W 
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Reexpressing the kinetic energy in terms of the momenta we find 
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kEB 
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The hamlltonlan is 


H = K + E 


where E la the configurational or potential energy discussed in section B 
The equations of motion for the atoms in the discrete region are 


R - 11 

5k, “ --.r 
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Thctse are just the standard classical equatlcns of motion. In addition, we 
alao find the new result ve are seeking, namely the equations of motion for 
the elastic field parweters. These are: 
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2H 
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where 

n ^ _L 0 li! p. 


and 


li . S’ v,e . 


1 , 5 ,,.. 


This specifies all the ingredients necessary to extend the kinetic energy 
quench method to the direct relaxation of the elastic field. It shows that 
each parameter of the elastic field obeys equations of motion similar to those 
for the atoms in the discrete region. In this sense these parameters behave 
as pseudo-atoms with inertia given by the mass tensor L , and subject to 
driving forces that arise partly from the potential energy, E, and partly 
from the elastic field dependence of the mass tensor. 

The presenr result shows that the ^'kinetic energy quench'' procedure 
may be applied to relax all regions of the model. It is a viable alternative 
to the use of the "conjugate gradients" approach of Fletcher and Reeves, In 


18 


fact there are scattered results based on the previous simulation approach 
indicating that the rate of convergence of the ''kinetic energy quench" pro- 
cedure nay be considerably faster than obtained from the Fletcher and Reeves 
approach. This probably is the case at least in the early stages of the 
relaxation, and suggests the sequential use of these two techniques for 
optimal convergence to equilibrium. 

Although the present extension of the equations of motions method 
la primarily designed as an aid to find the equilibrium positions, it may 
also contain the seeds for a bpoa fide modelling of time dependent crack, 
processes. 

E. The gradient of the configurational energy . 

The gradient of the energy turns out to be a sufficient and central 
ingredient in the two principal methods for relaxing the system to its equi- 
librium configurations (both to relative minlmums and to saddle points) . We 
show here the relationship of the gradient of the energy to the gradients of 
the individual pair potentials. 

The components of the gradient of the energy vector are the partial 
derivatives of the energy with respect to all the linearly independent para- 
meters in the energy. These are the positions of the atoms in the discrete 
region 

i, f. 3 


and the parameters specifying the elastic field which we shall denote by 


t ( 5 j • ■ • ; 


Three of these parameters we shall always take as the coordinates o£ the 
origin of the elastic field. Parameters beyond these correspond to non-linear 
eXatic field terms. 
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The gradient component 


In a discrete region atom coordinat e is 
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The gradient component In an elastic field parameter is related by 
the chain rule to gradients of the energy in the atom coordinates 
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The sum over discrete region atoms vanishes because of linear independence. 


i.e. , 
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The second sum ranges just over the atoms in C which interact with region D 
atoms, C*, because for all others in C we have 


- 0 


from the equilibrium condition assumed for the continuum. Therefore, we find 
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where the gradleats of the energy In a continuum region atom coordinates is 
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A particular comb Inat Ion of the gradients which is required in -the one-dimen- 
sional searches of the “ con j ugate gradient " method is the directional deriva- 
tive* This is 
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where the 


vector is an arbitrary vector, and is the magnitude thereof. 
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